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A   LATTICE   APPROACH   TO 
VOLUMES    IRRADIATED   BY 
UNKNOWN   SOURCES 

J.    Randa  and  M.   Kanda 

Electromagnetic  Fields   Division 
National  Bureau  of  Standards 
Boulder,    CO   80303 


We  suggest   an   approach  to  the   characterization  of   electromagnetic 
environments  irradiated  by  unknown  sources.     The  approach  is  based  on  the 
numerical  solution  of  Maxwell's  equations  subject  to  the  constraints   imposed 
by  the  measured  values  of  the  field  at   a  small   number  of  measurement   points 
and  by  boundary  conditions.      A  thorough  examination  of  two  methods  for   the 
numerical  solution   is   presented.      The  examples   attempted  demonstrate  the 
approach  but  reveal   that   neither   numerical   technique   is  fully  successful. 
Possible  future  directions   are  suggested. 


Key  words:      electromagnetic  environment   characterization;    electromagnetic 
environment  effects;    Hamilton's   action  principle;    ill-posed  problems; 
numerical  methods;    successive  over-relaxation  method. 


1 .    INTRODUCTION 

As  the  number  of  intentional  and  inadvertent  sources  of  electromagnetic 
(EM)  radiation  increases,  the  task  of  measuring  and  characterizing  EM 
environments  becomes  increasingly  difficult.   At  the  same  time  more  sensitive 
electronic  devices  are  being  used,  and  evidence  is  beginning  to  accumulate 
(see  e.g.  [1])  that  the  human  body  may  be  more  sensitive  to  low-level 
nonionizing  radiation  than  previously  assumed.   It  is  therefore  increasingly 
important  to  develop  methods  to  measure  and  characterize  the  EM  environments 
in  which  devices  and  people  are  found. 

The  only  systematic  approach  which  has  been  developed  and  actually 
implemented  is  the  statistical  approach  [2-4].   It  has  been  developed  quite 
thoroughly  for  applications  to  communications  systems  [5-7]  and  has  been  used 
in  a  wide  variety  of  applications,  including  EM  noise  in  mines  [8],  EM 
emissions  from  cars  [9,  10]  and  trains  [11],  radio  reception  inside  houses 
[12,  13],  and  characterization  of  urban  environments  [14-16].   Statistical 
methods  are  well  suited  to  cataloguing  general  features  of  general 
environments,  such  as  the  average  field  level  in  a  typical  room  of  a  generic 
suburban  house.   However,  they  require  too  many  measurements  to  be  efficient 
for  specific  characteristics  and  particular  situations.   (A  directional 
scanning  method  which  requires  less  measurement  effort  has  been  suggested 
[17],  but  it  has  not  yet  been  fully  developed.) 

In  this  paper  we  address  the  question  of  how  to  efficiently  extract 

information  about  a  specific  small  to  intermediate  (several  wavelength)  size 

environment,  e.g.  the  fields  in  one  particular  room  rather  than  the  average  of 

many  rooms.  We  assume  some  volume  of  interest  which  is  free  of  primary 

sources,  but  which  may  contain  conductors  carrying  induced  currents.   No 

knowledge  of  the  source (s)  irradiating  the  volume  is  assumed.   Instead  we 
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assume  knowledge  of  the  field  at  some  number  of  measurement  points.   In 

principle  the  field  measured  can  be  the  electric  or  magnetic  field,  or  both; 

but  in  this  paper  we  assume  it  is  the  electric  field.   If  both  were  measured, 

it  would  require  holding  fixed  spatial  derivatives  of  the  vector  potential  at 

measurement  points,  as  well  as  the  vector  potential  itself.  We  also  restrict 

our  attention  to  the  single-frequency  case.  The  extension  to  multiple 

frequencies  could  be  effected  through  superposition,  with  a  concomitant 

increase  in  measurement  effort. 

The  basic  idea  of  the  method  we  consider  is  to  attempt  to  numerically 

find  an  approximate  solution  to  Maxwell's  equations  which  is  consistent  with 

the  measured  values  of  the  field  at  measurement  points  and  which  is  also 

consistent  with  the  appropriate  boundary  conditions  at  the  surfaces  of  any 

conductors  within  or  bounding  the  volume  of  interest.   Although  it  is  not  the 

way  in  which  the  computation  proceeds,  one  can  think  of  the  problem  as  being 

divided  into  two  distinct  steps.   Since  the  geometry  of  a  real  situation  may 

well  be  quite  awkward,  the  normal  modes  of  the  system  would  not  be  known  and 

would  need  to  be  determined.  The  second  step  would  be  to  determine  what 

combination  of  modes  was  consistent  with  the  measurements.  This  second  point 

merits  some  discussion.  We  are  interested  in  obtaining  information  from  as 

few  measurements  as  possible,  and  therefore  in  the  typical  case  there  will  be 

insufficient  information  to  uniquely  determine  the  field  in  the  volume  of 

interest.   This  raises  two  complications.   Of  the  many  possible  solutions 

consistent  with  the  measured  values  of  the  field,  we  must  be  able  to  always 

find  the  "same"  one,  the  solution  with  some  identifying  feature  (preferably  a 

useful  one).   For  example,  we  would  want  to  always  find  the  one  with  the 

largest  (smallest)  average  energy  density,  or  the  most  probable  one,  etc.   The 

other  complication  is  that,  because  there  are  many  solutions,  the  matrix 
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equation  obtained  by  discretizing  the  system  involves  a  singular  matrix. 
Consequently,  methods  requiring  the  inversion  of  that  matrix  are  inapplicable. 

In  this  paper  we  develop  and  illustrate  the  approach  in  some  simple  cases. 
The  particular  geometry  chosen  is  a  rectangular  box  with  perfectly  conducting 
walls.   This  example  was  obviously  chosen  for  its  simplicity  rather  than  its 
practical  relevance,  but  it  is  not  so  far  removed  from  reality.   Introduction 
of  a  few  holes  for  windows  and  doors  makes  it  a  model  of  a  metal  building  or 
room.   Another  simplification  which  we  have  already  mentioned  is  the 
assumption  that  only  one  frequency  is  present. 

We  consider  two  different  numerical  techniques  for  finding  solutions  to 

Maxwell's  equations  subject  to  the  constraint  that  the  field  assume  specific 

values  at  the  measurement  points.   Both  methods  first  impose  a  spatial  grid  or 

lattice  on  the  volume  of  interest  and  are  formulated  in  terms  of  the  field 

values  at  the  lattice  points.   The  first  method  employs  Hamilton's  principle 

and  the  action  functional  for  classical  electromagnetism.   For  field 

configurations  which  are  solutions  of  Maxwell's  equations,  the  action  is 

stationary  with  respect  to  small  variations  of  the  field;  and  so  the 

appropriate  components  of  the  field  are  fixed  at  boundaries,  the  field  is 

fixed  at  the  measured  values  at  measurement  points,  and  it  is  allowed  to  vary 

everywhere  else  until  a  stationary  point  of  the  action  is  found.   The  second 

numerical  technique  is  to  attempt  a  direct  solution  of  the  relevant 

differential  equations,  which  are  the  vector  Helmholtz  equation  supplemented 

by  the  condition  that  the  divergence  of  the  vector  potential  vanish.   This  can 

be  considered  an  alternate  approach  to  locating  the  stationary  point  of  the 

action,  since  the  action  stationarity  leads  to  Maxwell's  equations  and  thence 

to  the  Helmholtz  equation  and  divergence  condition  in  our  case.   Again  the 

field  is  set  equal  to  measured  values  at  measurement  points,  and  the  Helmholtz 
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equation  and  divergence  condition  are  solved  simultaneously  using  a  successive 
over-relaxation  (SOR)  method.   Neither  of  the  two  numerical  methods  is 
entirely  successful,  but  the  SOR  method  has  enough  success  to  demonstrate  the 
general  approach  and  to  point  out  avenues  for  development. 

In  the  next  section  we  present  the  first  numerical  technique,  based  on 
the  action  functional.   The  SOR  technique  is  applied  to  the  Helmholtz  equation 
and  divergence  condition  in  Section  3.   In  the  final  section  we  summarize  the 
work  and  comment  on  future  prospects. 

2.    STATIONARITY  OF  THE  ACTION 
2.1    Formulation 

In  terms  of  the  usual  vector  and  scalar  potentials  A(x*,t) 
and  cj>(x,t),  the  action  functional  for  classical  (lossless)  electromagnet  ism  is 
[18,  19] 

,t2 


&=   Jt  dt  Jd3x  {![(%  +  f^  A)  •   e(x)  -  (vcj,  +  §£  A) 

-  (V  x  A)  •  y_1(x)  •  (V  x  A)  +  A  •  J(x,t)  -  <j>  p(x,t)}, 

where  "?  and  ti  are  the  permittivity  and  permeability,  and 

J  and  p  are  the  current  and  charge  densities.   The  spatial  integral  extends 
over  the  volume  of  interest,  subject  to  the  restrictions  imposed  by  the 
variations  of  the  potentials  vanishing  on  the  surface  (see  below).   We  work 
with  the  potentials  rather  than  the  electric  and  magnetic  fields  because 
Hamilton's  principle  requires  that  the  quantities  being  varied  be  independent 
degrees  of  freedom,  whereas  E  and  H  are  related.   The  four  potentials 
(A,(}))  cannot  all  be  independent  either,  of  course,  but  one  can  be  easily 
eliminated  when  the  gauge  choice  is  made. 
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In  this  formulation,  the  two  homogeneous  Maxwell's  equations  follow 


-> 


immediately  from  the  definition  of  E  and  B  in  terms  of  the  potentials: 

E  =  -  (v^  +  f-  a)  I  ^xE+|-B  =  0 

3t      ►  3t  (2.2) 

B  =  ^xA  '  V  •  B  =  0 

The  two  remaining  Maxwell's  equations  follow  from  Hamilton's  principle  [10]. 
If  we  consider  small   variations  of  the  potentials  and  requires  that  the 
resulting  variation  of  the  action  vanish,  we  obtain 


(2.3) 


6ACl  =   0   *     V   •   D  =   p, 

6    n    =   0   +     7  x  H  -  |-  D  =  J, 

A  8t 

-> 

provided  the   variations   of  A   and   <j>  vanish  at   the   end  points.      To  be  more 

exact,   what  is   required   is   that   6<j>  and  6A   vanish  on  the  surface  of  the   volume 
V  at   all   times   and  that   6k  vanish  throughout   the   volume  at   times   tx   and  t2. 
What   this  means   in  practice   is   that   in  order   to  use  this  in  a  calculation,   we 
must  specify  the   fields  on  the   boundary  surface  at   all   times  and  throughout 
the   volume   at   the  initial   and  final    times.      In  any  real   application,   such  a 
superabundance  of   information  will   not    be   available  for   the   volume  of   interest, 
In  order  to  exploit   the  stationarity  of  the  action  we  must  expand  the  volume 
considered  beyond  the  region  of   interest,   out   to  distances  where  the   fields 
can   be   assumed  to   be  negligible.      The  same  applies   to  the   time;    t1    is   chosen 
before  the   fields   are  turned  on  and  t2   after   they  are  turned  off. 


For  the  calculations  that  follow,  it  will  be  convenient  to  work  in  a 
definite  gauge.   A  propitious  choice  is  cj>(x*,t)  =  0.   The  action  then 
assumes  the  form 

d-   Jdt  jd3x  |l[(f-  t)   •  %   •  (|r  t)  -   (?  x  1)   .  ?-«  •  (v-  x  1)] 

2   9t  3t  (2>4) 

+  A  •  J}. 

This  gauge  has  the  convenient  feature  that  for  the  single-frequency  case  a 
measurement  of  the  electric  field  directly  determines  A, 
E(x,t)  =  a)  A  (x,t  -  ir/2w). 

Having  extended  the  volume  under  consideration  to  virtually  all  of  space- 
time,  we  must  restrict  the  problem  to  manageable  size.   In  particular,  since 
we  do  not  know  the  sources,  the  volume  of  interest  is  chosen  to  exclude  them. 
We  envision  the  division  depicted  in  Fig.  2-1.  The  volume  marked  Vj  is  the 
region  of  interest;  it  is  assumed  to  be  free  of  primary  sources,  but  it  may 
contain  conductors  with  induced  currents.  Volume  V3  is  a  buffer  zone 
separating  Vj  from  V5,  which  is  the  rest  of  space,  wherein  are  located  any 
primary  sources.   The  idea  is  to  divide  the  action  into  one  piece  from  the 
integral  over  volumes  Vj  and  Vg  and  another  piece  from  V5.   Because  Vg 
contains  unknown  sources,  we  do  not  try  to  determine  the  fields  there,  which 
means  that  we  also  will  not  know  the  fields  on  the  surface  between  V5  and  Vg. 
The  fields  on  that  surface  will  be  allowed  to  vary  or  will  be  fixed  by  a 
reasonable  guess.   Obviously,  near  the  boundary  between  Vg  and  Vg  the 
solution  obtained  will  be  very  sensitive  to  the  choice  for  the  fields  on  that 
boundary,  and  therefore  it  will  not  be  reliable.   As  we  move  away  from  the 
outer  boundary  of  Vg  the  values  of  the  fields  should  be  influenced  more  by  the 
measurement  points  and  less  by  the  values  on  the  surface  between  Vg  and  Vg. 

For  points  far  enough  away  from  that  surface — i.e.  within  VT  —  it  should  be 
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possible  to  obtain  reliable  solutions,    given  enough  measurement   points.      The 
hope   is   that   "far   enough  away"   and   "enough  points"   are  not  so  large   as   to 
render  the  method  impractical  for  most  applications.     The   positions  of  the 
measurement   points  will   clearly  affect  the  size  required  for  Vg;    it  may  well 
be  advantageous   to  make  a   few  measurements   on  the   perimeter   of  Vg.      It  would 
probably  also  be  advisable  to  choose  the  boundary  of  Vg  to  coincide  with 
conducting  walls  when   it   is  feasible,    in  order   to   constrain  the   fields  on  the 
boundary  as  much  as   possible. 

The   quantity  we  shall   consider  then   is  a  reduced  action,  (Z,   which   is 
defined  as  in  eq.    (2.4)   but  with  the  integral  restricted  to  Vj  €>  Vg.      (In  the 
general  multiple-frequency  case  the  volumes  are  four-dimensional  space-time 
volumes.)     The  next  step  is   to  discretize  the  expression  for  the  reduced 
action,    converting  the   volume   integral   into  a  summation  which  approximates   it. 
There   are   any  number   of  ways   to   do  so;    we   are  not  interested  in  their  relative 
(dis)advantages  at   this  time.      The  discretized  reduced  action  will   have  the 
general  form 


y\ 


a 


dt 


I 
a,B,YeVI@V{ 


AV 


1 


3      * 


aeriKft^'U-c, 


BY 


(ft*(t)U 


(2.5) 


-   tf.   t(t»  •   M-Y   •    (*"  tCt»        I    • 


The  discrete   indices   a,B,Y  label   the  spatial   points,   on  which  are  centered  the 
volume  elements   AVagy      We  have   left   the   time  variable  continuous  for   now, 
anticipating  the   single-frequency  example  below.      The   quantities 
(3A/8t)agy  and    (V   *   A)agY  will    in  general   depend  on  the   values  of  the 


field  A  at  a  number  of  points  on  (or  within)  the  surface  bounding  AVagy 
A  concrete  realization  of  this  discretization  of  &  will  be  given  in  the 
example  below.   Imposition  of  the  constraints  required  by  the  known  values  of 

E  and/or  H  at  measurement  points  can  be  rather  complicated,  but  for  only  one 

->      -> 
frequency  and  for  our  gauge  choice  Eagy  <*  Aagy,  and  a  measurement  of 

->■ 
the  electric  field  at  a  point  fixes  Aagy  at  that  point.   A  similar  comment 

applies  to  boundary  conditions  at  perfectly  conducting  walls;  for  the 

single-frequency  case  they  can  be  imposed  with  relative  ease. 

The  calculation  then  proceeds  as  follows.   A  grid  is  defined  within 

Vg@Vj,  and  Aagy  is  fixed  at  measurement  points  and  appropriate 

-> 
components  are  set  equal  to  zero  at  conducting  walls.   All  other  A^y's 

are  varied,  and  we  search  numerically  for  a  stationary  point  of  (X   in 

eq.(2.5).   When  (if)  a  stationary  point  is  found,  then  that  set  of  Aagy 

constitutes  an  approximate  solution  to  Maxwell's  equations  within  VggVj 

which  is  consistent  with  measured  field  values  and  boundary  conditions.  The 

problem  of  finding  the  stationary  point  is  nontrivial,  however,  and  we  now 

turn  our  attention  to  that. 

2.2    Nature  of  the  Stationary  Point  and  the  Euclidean  Action 

The  technique  presented  in  this  paper  depends  on  finding  field 

configurations  such  that  the  action  is  stationary  with  respect  to  small 

variations  of  the  field.   Such  field  configurations  are  solutions  of  Maxwell's 

equations.   The  standard  method  of  finding  the  stationary  point  would  be  to 

discretize  the  functional  as  in  (2.5),  set  3#./3A^gy  =  0  to  obtain  a 

(large)  matrix  equation,  and  solve  that  equation.   That  approach  founders  here 

because  the  matrix  would  be  (very)  singular.   There  are  in  general  many 

possible  solutions  consistent  with  boundary  conditions,  measured  fields,  and 

unknown  sources.  We  need  to  locate  a  stationary  point  and  to  know  which 
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stationary  point  it  is — does  it  correspond  to  the  solution  with  the  largest 
maximum  electric  field  in  the  volume,  or  the  solution  with  the  minimum  power 
density,  etc.   To  locate  the  stationary  points  numerically  we  need  to  know 
whether  they  are  maxima,  minima,  or  neither.   The  phrase  "Principle  of  Least 
Action"  may  lead  one  to  expect  a  minimum,  but  in  general  this  need  be  true 
only  for  infinitesimal  time  intervals  At  =  t2  -  tx.  To  consider  our  specific 
case,  we  adapt  the  treatment  of  Whittaker  [20].   We  use  the  form  of  eq.  (2.4) 
and  assume  for  simplicity  that  *£   and  "ft-1  are  scalar  constants.   If 
the  fields  are  allowed  to  vary 


A  (x,t)  ■»  A(x,t)  +  ct(x,t), 


a(x  ,t  j )    =   a(x  ,t2)    =   0, 
where  3  is  small,   then  the  second-order  variation  of  the  action  is 


(2.6) 


62a  =  \  , 


t2 


dt 
t, 


d3x    [ea2   -  -  ($  x   a)2].  (2.7) 

V 


The   usual   argument   is   that   because  5  but  not   3  vanishes   at   tx    and 

t2,    for  small   t2-tx    the  first   term   dominates  the   integral   in    (2.7)    and  the 

stationary  point   is   a  minimum.      That   argument   fails   in  a   field  theory, 

however,   where  the  second  term  involves  spatial   derivatives  of  the  dynamical 

variables.      Because  the   functional   form  of  &(x*,t)   is  arbitrary,    its 

spatial   dependence  can  oscillate  wildly  and  62CL  can  be  positive  or  negative 

(depending  on  the   choice  of  3)    for   any  fixed  nonzero   At.      Therefore  the 

stationary  points   of     are  neither  maxima  nor  minima. 

This   changes  when  the  action  is  discretized  by  the  introduction  of  a 

spatial   grid.      Then  the  spatial   derivatives  of  dt  are   bounded  by 

(V   x   &) 2  <   0(&2/a2),    where   a   is   the   lattice  spacing.      The  stationary  point 

is   a  minimum  provided   At   is  less   than  some  number  on  the  order  of  /e\Pa.      For 
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larger   values  of  At,    62d  can  again  have  either  sign,   and  consequently  the 
minimum  is  again  a  saddle  point.     That  is  the  case  of  interest  since  we  need 
At  ■+  °°  in  order  to  set  all   the  fields  equal   to  zero  at   tx   and  t2. 

This  poses  a  serious   calculational  problem.      A  stationary  point  which  is 
neither   a  maximum  nor  a  minimum  will   be  very  difficult  to  locate  numerically. 
In  earlier  work   [21,    22]  we  introduced  a  method  which  succeeded  in  some  simple 
cases,   but  subsequent   tests  with  more  complex  field  configurations    (i.e. 
higher  modes)   have  exposed  inadequacies  in  that  method.      Faced  with  the 
impracticality  of  numerically  locating  saddle  points  of  functions  of  very  many 
(103-106)   variables,   we  need  a  function  to  maximize  or  minimize.      Our  efforts 
to  construct  such  a  function  have  not   been  marked  by  complete  success.      One 
interesting  attempt  which  was   partially  successful  is  to  borrow  a  trick  from 
quantum  field  theory  and  work  in  Euclidean  space-time   [23].      To  do  so  we  take 
time  to  be  imaginary,    t  =  -  it  with  t  real,   thereby  changing  the  sign  of  the 
(9A/8t)2  term  in   (2.5)   and   (2.7).     We  solve  the  problem  in  Euclidean  space  and 
analytically  continues  the  result   back  to  Minkowski  space,    x+it.      The 
Euclidean  action  is 


V 


d3xfl(|_XE(;,T))  .t.(f^E) 


dx 
Tl  (2.8) 


+  \   (v-   x   AE)    •    y"1    •    (v-   x   AE)    -  AE    •    J}    , 

and  this  is  to  be  minimized  to  obtain  the  physical  solutions  in  Euclidean 
space.  The  analytic  continuation  back  to  Minkowski  space  can  be  a  difficulty 
in  general,  and  in  fact  it  will  prove  to  be  the  limiting  factor  for  this 
technique.  Nevertheless,  it  is  possible  to  handle  some  simple  examples  in 
this  manner. 


11 


As  noted  above,  there  will  be  many  possible  solutions  consistent  with  a 
few  field  measurements  if  the  sources  are  unknown.  The  Euclidean  solution  is 
obtained  by  minimizing  the  integral  of  the  energy  density 
(V2  e  Eg  +  V2  \x   Hg).   As  will  be  seen  below,  for  the  single  frequency  case 
this  also  corresponds  to  minimizing  the  energy  density  in  Minkowski  space. 
Therefore,  the  solution  obtained  will  be  the  one  with  the  minimum  average 
energy  density  in  the  volume  under  consideration.  This  is  a  welcome  bonus  of 
working  in  Euclidean  space--not  only  can  we  find  a  solution,  but  it  is  a 
useful  one  as  well. 

The  entire  procedure,  including  the  analytic  continuation  and  numerical 
search  method,  is  greatly  clarified  by  an  example,  to  which  the  next  sub- 
section is  devoted. 

2.3    Simple  Example 

A.   Problem  and  Calculation 

Having  presented  the  general  ideas  of  this  approach,  we  now  attempt  to 

implement  it  in  a  simple  example — a  rectangular  waveguide  with  perfectly 

conducting  walls.   This  is  obviously  not  supposed  to  be  a  practical 

application,  and  many  of  the  difficulties  and  nuances  of  the  general  case  are 

absent.   It  is  a  practice  problem  to  demonstrate  the  idea  and  to  provide  a 

basis  on  which  to  build  toward  solution  of  real  problems.  The  rectangular 

waveguide  is  chosen  because  it  has  simple  known  modes,  because  the  boundary 

conditions  are  easily  imposed,  and  because  it  can  be  reduced  to  a  two-(or  even 

one-)  dimensional  problem,  thereby  reducing  the  computational  exercise.  As  we 

shall  see  below,  even  in  this  simple  case,  it  is  not  entirely  trivial  to 

obtain  the  correct  continuum  answer,  with  which  to  compare  the  results  of  our 

numerical  computation. 
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We  first  restrict  the  problem  to  two  dimensions,  obtaining  the  general 
forms  for  Ag  and  A^  (the  vector  potentials  in  Euclidean  and  Minkowski  space, 
respectively)  and  the  relationship  between  them.  To  facilitate  the  analytic 
continuation,  we  assume  a  standing  wave  due  to  perfect  reflection  from  the 
(very  distant)  end  of  the  waveguide.   For  a  given  TE  mode,  TEmn,  A^  can  then 
be  written 

A  (x,t)  =  A(xi)  cos  kz  cos  cot  =  A(]_)  cos  kz  cos  tot, 


(2.9) 


/to2    2  rm2   nS 


where  a  and  b  are  the  transverse  dimensions  of  the  waveguide    (aSb),    c   is  the 
speed  of  light  in  the  waveguide  medium,   xl    =   (x,y)    is  the  two-dimensional 
transverse  position  vector,    and  A^(x,t)    and  A(J_)    are  real.      Equation   (2.9) 
would  allow  a  simple  restriction  to  two  dimensions  for  a  single  mode. 
However,   the  solutions  of  interest  are  those  with  the  minimum  energy  density 
consistent  with  the  measurements,   and  those  solutions  are  superpositions  of 
many  modes.     Consequently  the  true  solution  will   not  have  a  z  dependence  which 
can  be  factored  out  as  in   (2.9),   and  the  problem  will   not   reduce  to  two 
dimensions.      In  order  to  force  it  to  do  so,   we  impose  the  form  of  eq   (2.9) 
with  k  taken  to  be  that  of  the   lowest  mode  consistent  with  the  measurements. 
That  yields  a  two-dimensional   problem,   but  it  is  no  longer  electromagnetism: 
if   the  solutions   A(J_)    are  substituted  back   into    (2.9),    the  resulting 
A(x,t)   need  not   be  a  solution  of  Maxwell's  equations.      Nonetheless,    the   action 
and  the  problem  are   quite  similar  to  electromagnetism,   and  it   is  a  suitable 
problem  for  demonstration  purposes  since  one  can  obtain  the  true  solutions  for 
comparison.     This  arbitrarily  imposed   z  dependence  does  not  reflect  a 
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shortcoming  of   the  method,   but  rather   the   difficulty  of  finding  a  simple 
two-dimensional   application. 

To  obtain  the  general  form  for  Ag,   we  require  that   it  be  real  and  that 
it  satisfy  the   Euclidean  wave   equation 

(V2    +   Ue  |^r)   Ae(x,t)    =   0    .  (2.10) 

It  can  then  be  written  as 


A_   (x,t)    =  aR(_L)    e  +  b„(J_)    e 


(2.11) 


■>   ,iv      -i(kz   +  ojt)       *   ,h      -i(kz  -  wx)  . 

+  cE(J_)    e  +  dE(J_)    e  +  compl.    conj . , 

where  u  need  not  be  (in  fact  is  not)  real.   If  (2.11)  is  then  continued 
back  to  Minkowski  space  (x  =  it)  and  compared  to  (2.9),  the  relation  between 
Ag  and  A^  follows.  We  can  write 

A   (x,x)  =  A(J_)  cos  kz  cosh  ux ,  (2.12) 

where  A(J_) ,  k,  and  oj  are  the  same  as  in  (2.9).   Equation  (2.12)  enables  us  to 
use  measurements  of  A^  to  constrain  Ag.   Conversely  once  Ag  is  determined  by 
the  computation,  A(J_)  and  hence  also  A^  are  known. 

Substituting  the  form  for  Ag  (2.12)  into  the  Euclidean  action  (2.8)  leads 
to 

a+A  b  +  A 


aE  =  c 


dx 

-A  -A 


U)2 


dy  Kpr  +  k2)  [a2(1)  +  a2(D]  +^a2(1) 


(2.12) 


+    [(9  A    (I))2    +    (9  A    (I))2    +    (3   A    (]_)    -    3  A    (I)]2] 
y   z  -1-   '  x   z  -1-   '  v    x  y  -1-  y  x  -1-   ;    J 

-2y   J(l)     •    A(]_)}, 
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where  J(j_)    is   defined  in  terms  of  the  wall  current  J(x,t)    in  the 
obvious  way   (9),   and  where  the   constant   C  is   given  by 


T 


2 


2U    . 


dx  cosh  ut 


dz  cos2kz.  (2.14) 


In  writing  (2.13)  we  have  assumed  that   the  material  in  the  waveguide   is 
isotropic  and  that  the  range  of  the  z  integration  is  either   very  long  or   an 
integral  number   of   cycles.      The   induced  current  and  the  limits  of  the 
transverse  integrations  require  explanation.      In  applying  Hamilton's   principle 

the  variations  of  A  must   be  zero  on  the  boundary,  which  in  our  two-dimensional 

->■ 
case  here  means   that  we  must  specify  A  on  the  transverse  boundary.      In  order 

-> 
to  be  able  to  specify  all   components  of  A  we  choose  the  boundary  to  lie  a  few 

skin  depths  within  the  conducting  walls  where  the  field  can  safely  be  assumed 

to  vanish.      In  figure  2-2  the  x  integration  goes  from  x  =   -A  to  x   =  a  +   A,    and 

the  y  integration  from  -A  to  b  +  A,   where  A  =  N6,    some  suitable  number  of  skin 

depths.     Then,   however,   the  currents   induced  in  the  walls  are  contained  within 

the   integration  volume.      Fortunately,  we  can  show  that   the  contribution  to  the 

action  from  the  volume  within  the  conductor  is  negligible,   and  we  can  write 

a  b 


aE  =  c 


dx 


dy  !(£  +  k2)  [a*<]_)  *  a2(1)]  +^(1) 


(2.15) 


+    [(3      A    (I))2    +    (3  A    (I))2    +    (3  A    (I)    -    3  A    (I))2]}. 

L^yzJ_;  (.  x  z  J_  J  VxyJ_  yXi_;JI 

We  still   need  to  impose  Etan  =   0  and   Bnorm  =   0  at   the   conductor  walls,   which 

is  accomplished  by  the  requirements  that  A(]_)^an  var>ish  at   the  walls. 

The  action  is  then  discretized  by  introducing  a  rectangular   grid  into  the 

waveguide   as   indicated   in  Fig.    2-3.      The  spacing  between   points   is   Ax   =  a/Nx, 

Ay  =  b/Ny.     Derivatives  at  a  point    (i,j)   are  defined  by 
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[j_A)ij.ian..i)-A(i„i)t 

except  at  the  boundaries,   where 

X  X  , 

and  similarly  for   y  derivatives.      The  Euclidean  action  then  takes   the  form 

N        N 
/7  x       y 

"-„     =  C  X        I       AA..    {[A   (i,j)2   +  A   (i,j)2]   c     tt2 
E  .    .    .    _        lj    lL  x     'J  y     '°      J     mn 

i=0  j=0 

+   4it2   A   (i.J>«   +    [(A  (l.J  +  1)    -  A   (i,j))/Ay]2 

Z  Z  Z  (2.18) 

+    [(Az(I+l,j)    -  Az(i,j))/Ax]2    +    [(Ax(i,j  +  1)    -  Ax(i,j))/Ay 
-    (Ay(i+l,j)    -  Ay(i,j))/   Ax]2}, 

Cmn  =    ^8   "  f*  "  ^    ' 

where  all  lengths  are  in  units  of  the   free-space  wavelength  A.      The  area 
elements   are 

AA        =   Ax    Ay  0<k,KN 

=  V-jAx   Ay  k  or   1   =   0  or  N,  (2.19) 

and  the   "out  of  bounds"   fields  are  defined  as 


(2.20) 


A(i,Ny   +   1)    =   2A(i,Ny)    -   A(i,Ny  -    1), 

A(NX    +    1,    j)    =   2A(Nx,j)    -   A(NX   -    1,j). 

Having  discretized  the   problem,   we  next  set   the   values   of   the  tangential 
fields   to   zero  on  the   boundaries,    set   the  fields  equal   to  their  measured 
values   at  measurement    points,    and  vary  all   other   fields   to  minimize  6cg   in 
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(2.18).      The  manner   in  which  we  minimize^   is  to  step  through  the  grid,   at 
each  point    (i,j)    setting  the  fields  equal   to   the   value  required  by 
8  <3.g/3Aa(i  ,j )    =  0,    a   =   1-3,    given   the   current   values  of   the   field  at 
neighboring  points.     This   process  is  repeated  until  the  action  reaches  a 
minimum  and  does  not   change  significantly  during  further   passes   through  the 
grid.      (This  is  just   the   Gauss-Seidel  method;    see  e.g.    [24].) 

B.   Demonstration 

To  further   demonstrate  the  method  and  to  clarify  the  sort  of  solution 
obtained,   we  present  a  simple  one-dimensional   computation  and   compare  it   to 
the  continuum  results.      The  field  was  fixed  at   three  measurement   points, 

Ay    (|,|)    =  0.707,    Ay    (|,|)    =   1.0,    Ay    (-jp,|)    =   0.707,  (2.21) 

with  Ax  and  Az  equal  to  zero  at  all  three  points.  This  field  pattern 
corresponds  to  the  TE10  mode,  or  any  number  of  higher  modes.  The  waveguide 
dimensions  were  taken  to  be  a  =  2. 0,  b  =  1 . 0.   The  computation  to  find  the 
minimum  is  straightforward  and  proceeds  as  outlined  at  the  end  of  the  previous 
subsection.  Obtaining  an  analytical  result  for  comparison  is  less  trivial. 
To  do  so  for  the  three-measurement  case,  we  assume  the  y  dependence  of  all 
fields  is  a  constant  and  expand 

00 

Ay(x)  =  I     a.  sin(j~).  (2.22) 

J-1   J 

Ax  and  Az  are  zero  at  the  measurement  points  and  are  therefore  zero  everywhere 
for  the  minimum- UE  solution.   In  terms  of  the  aj  of  (2.22)  the  action  is 

abTT2 

VCT~      Ia?31    +^2)-  (2.23) 

J      J 

The  constraints  of  (2.21)  can  be  combined  to  yield 
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ai  +   I   <a8n+1  "  a8n-l)  =  1» 
n=1 


(2.24) 


and  two  other  constraints  which  lead  to  all  a^  which  do  not  contribute  to 
(2.24)  being  zero. 

In  order  to  minimize &e   (2.23)  subject  to  the  constraint  (2.24)  we  first 
use  the  constraint  to  eliminate  ax  and  then  set  derivatives  of  (Xe   with  respect 
to  all  other  a^  equal  to  zero.   This  yields  the  infinite  dimensional  matrix 
equation 


AT 

1 
1 
1 


At 


1 

1 

a; 
i 


1  . 
1  , 

1  . 

At 


bT 


,(2.25) 


where 


bn  =  ±  a8n±1 


A^  =  2  +  2n2  ±  n/2, 


(2.26) 


This  system  can  be  solved  by  subtracting  each  row  from  the  first  to  get 
the  relation 


±    A,  -  1 

)   =  —i 

An  ~  1 


bx. 


(2.27) 


The  first  row  of  (2.25)  then  determines  b7,  and  ax  can  be  determined  from 
(2.24).   The  result  is 


a8n±1  =  ±  — 


1 


ax ,  n>0, 


(2.28) 


An  "  1 
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..  -  [1  ♦    I  (-^-♦^rJ-)] 


- 1 


n=1  An  -  1    An  ■ 


4tt 


=  j^   tahn  [1  /3T1)  =  0.43530, 


Equation  (2.28)  and  the  fact  that  all  other  a^  =  0  can  then  be  used  in  (2.22) 
to  determine  the  true  solution  to  which  to  compare  the  result  of  the  numerical 
computation.   Figure  2-4  shows  the  comparison,  and  it  is  obvious  that  a)  the 
two  sets  of  results  agree  with  each  other,  and  b)  they  agree  with  the  measured 
values  at  x  =  a/4,  a/2,  and  3a/4. 

2.4  Comments 

The  simple  application  in  the  preceding  subsection  illustrates   both  the 
general  lattice  approach  to  complex  environment  characterization  and  the 
attempt  to  implement  this  approach  using  Hamilton's  principle  and  the  action 
in  Euclidean  space-time.      The  example  demonstrates  that  the  lattice 
computation  does   produce  the  correct  results.      It  indicates  that   in  one 
dimension,   for  a  system  very  similar  to  Maxwell's  equations,   we  are  able  to 
find  the  minimum-energy-density  solution  consistent  with  measured  field  values 
at   a  few  points,  without   knowledge  of  the  external  sources. 

Unfortunately,  a  fundamental  obstacle  arises  when  less  trivial  examples 
are  attempted.  If  there  is  any  spatial  variation  in  the  phase  of  the  field, 
or   if  the  different   components  of  the  vector  potential   are  not  all   in   phase 

with  each  other,    it  becomes   impossible  to  perform  the  analytic  continuation  as 

->■ 
required.      Since  we  are  taking  A   to   be  real,    phase  variation  appears   as   an 

admixture  of   cos   mt  and  sin   wt   time   dependences.      One   cannot   continue   both   cos 
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and  sin  to   imaginary   time   and  still  maintain  the  reality  of  A   (which  is 
necessary  in  order   to   keep  the   Euclidean   action  positive   definite). 
Consequently,    it  does  not  appear  that  the  trick  of  analytically  continuing  to 
imaginary  time  will   be  applicable  in  practical   problems. 

That  leaves  us  with  the  problem  of  finding  a  stationary  point  of  the 
action  which   is  neither  a  maximum  nor  a  minimum.     Until   that  numerical  problem 
can  be  solved,   the  action-based  approach  holds  little  promise. 

3.         DIFFERENTIAL  EQUATION   SOLUTION 
3.1  Method 

As   in  the   preceding  section  we  assume  a  single  frequency,    and  we  continue 
to  choose  the  gauge  so  that   the  scalar   potential  vanishes.     We  then  have 

<j>(x,t)    =  0,        A(x,t)      =     A(x)    e"lut, 

E(x)    =   ia)A(x),        kx)      =     V   x   fox).  (3'1) 

Equation  (3.1)  guarantees  that  E  and  B  will  satisfy  the  two  homogeneous 
Maxwell  equations.   In  a  sourceless  region  of  free  space  the  two  remaining 
Maxwell  equations  can  be  written 

V  •  A(x)  =  0,  (3.2a) 

V2  A(x)  +  o)2ue  A(x)  =  0.  (3-2b) 

We  continue  to  measure  all    distances   in  units   of   the   free  space  wavelength 
A   =   2tt/(u>  /pe),    so  that   w2ye   =   4tt2    in    (3-2b). 

The  method  presented   in  this   section   is   to  discretize  the  differential 
equations   of    (3.2)    and  then  attempt   to  numerically  solve  them  simultaneously 
subject   to  the   constraint   that   the   field   is  fixed  at  measurement   points. 
Using  finite   difference   forms  for   the   derivatives,    eqs    (3.2a)   and    (3.2b) 
become 
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Ax  ( i  + 1  ,  j  ,  k )    -Ax(i-1,j,k)    +Ay(i,j+1,k)    -Ay(i,j-1,k)    +Az(i,j,k  +  1) 

-Az(i,j,k-1)    =  0  (3.3a) 

[A(i+1,j,k)    +  A(i-1,j,k)    +  A(i,j+1,k)    +  A(i,j-1,k)    +  A(i,j,k  +  1) 

+   A(i,j,k-1)    -   6    A(i  ,j  ,k)  ]/A2    +    4rr2    A(i,j,k)    =   0,  (3. 3b) 

where  we  have   introduced  a  rectangular   grid  with  equal   spacing   A   in  x,y,    and 
z,   and  the  positions  of  grid  points  are  given  by    (x,y,z)    =   (i,j,k)A.      As  in 
the   previous  section  the   cavity  will   have   dimensions   a,b,c,   and    (Nx»Ny,Nz)    = 
(a,b,c)/A  are  the  number   of   intervals   in  each  direction.      Points   on  the 
boundaries   require  special   treatment   since   they  lack  neighbors   in  one  or  more 
directions.      There  we  use  the  equation   at   the  neighboring   interior   point.      For 
example,    the   field  at   a  point    (a,y,z)    =   (Nx,j,k)A   is   determined  from 
Ax(Nx,j,k)    -   Ax(Nx-2,j,k)    +  Ay(Nx-1,j+1,k)    -   Ay (Nx-1 , j -1 , k) 

+  Az(Nx-1,j,k+1)    -   Az(Nx-1,j,k-1)    =   0,  (3.4a) 


[A(Nx,j,k)    +  A(Nx-2,j,k)    +  A(Nx-1,j+1,k)    +  A(NX-1 , j-1 ,k)    +  A(NX-1 , j , k+1 ) 
+   A(Nx-1,j,k-1)    -   6A(Nx-1,j,k)]/A2    +   4tt2   A(Nx-1,j,k)    =   0.      (3.4b) 


In  general  there  will  be  more  than  one  possible  solution.   A  unique 
solution  would  require  measurements  at  sampling  theorem  separations  (<  A/2  for 
infinite  free  space),  or  knowledge  of  the  sources,  or  knowledge  of 
A  on  the  boundary  surface  of  the  volume.   This  restricts  the  range  of 
numerical  methods  available.   To  see  this,  consider  just  the  vector  Helmholtz 
equation,  (3.3b).   By  arranging  all  field  components  Ax(i,j,k),  Ay(i,j,k), 
Az(i,j,k)  into  one  long  vector  z,  one  can  put  (3.3b)  into  the  form  of  a  matrix 
equation  Bz  =  0.   Those  elements  of  z  corresponding  to  measurement  or  boundary 
points  are  known  and  can  be  brought  over  to  the  right  hand  side  to  yield  an 
equation  of  the  form  Ax  =  b.   If  there  is  more  than  one  solution,  then  A  is 
singular  and  any  method  requiring  the  existence  of  A-1  will  fail. 
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This  leads   us   to   adopt  an   iterative  technique — the  successive 
over-relaxation    (SOR)  method   [24].      For   simplicity,    consider  first  just   the 
vector  Helmholtz  equation,    (3. 3b).      To  apply  SOR   to  this  we  step  through  the 
grid,    and  at   each  point    (i,j,k)   we   first   compute  the  residual,    defined  by 

RH(i,j,k)    =  A(i,j,k)    -    [A(i+1,j,k)    +  A(i-1,j,k)    +  A(i,j  +  1,k)  ,        . 

■>  *  +  (3'5) 

+  A(i,j-1,k)    +  A(i,j,k  +  1)    +  A(i,j,k-1)]/(6   -  4u2A2). 

(This   is  actually  the  residual   divided  by    (6   -  4ir2A2),    but  we  shall   call   it 
the  residual.)     We  then  change  the   vector   potential   at   that   point  according  to 

A'(i,j,k)    =  A(i,j,k)   -   ftH  RH(i,j,k),  (3.6) 

where   &#   is   a  number   between   0  and   2,    chosen  to  optimize  the  rate  of 

convergence.      For  Q,^  =   1  ,    this   is  the  Gauss-Seidel  method,    in  which  the 

-> 
function    (A)    is   set   equal   to  the   value  required  to  satisfy  the   difference 

equation(3«3b) ,    given  the  current   values  of  the  function   at   neighboring  points, 

For   fiH  greater    (less)    than  one,    the  method   is  one  of  over-(under-)   relaxation. 

The  term  SOR   is   used   generically  to  apply  to   both  over   and   under   relaxation 

(and  Gauss-Seidel   as  well)    [24], 

Application  of   the   SOR   method  to  the   divergence   condition    (3. 3a)    is 
somewhat   different   due  to  the  fact   that    (3.3a)    relates   different   components  of 
A  at   different    points,    and  we  would  like  to  treat   all   three  components   the 
same.      This   can   be  accomplished   in  the  following  manner.     We   define  the 
residual   as 

Rv(i,j,k)    =   -Ax(+)    +  Ax(-)    -   Ay(+)    +  Ay(-)    -   Az(+)    +  Az(-),  (3.7) 

where 

Ax(±)    =  Ax(i    ±  1 ,j  ,k), 

Ay(±)    =  Ay(i,j    ±  1,k),  (3.8) 

Az(±)    =   Az(i  ,j  ,k    ±   1 )  . 
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The  relaxation   is  then   performed  by   spreading  the  residual   equally  over   the 
Aa(±)'s, 

A^(±)    =   Ax(±)    ±  \     fty  Rv(i,j,k), 

A*(±)    =  A   (±)    ±  1     n7Rv(i,j,k),  (3.9) 

A^,(±)    =  Az(±)    ±  \     Qv  Rv(i,j,k), 

where  the  relaxation  parameter  fty  e(0,2)  is  again  chosen  to  optimize 
convergence. 

We  have  discussed  separately  the  SOR  method  for  the  vector  Helmholtz 
equation  and  for  the  divergence  condition,  but  in  practice  we  want  to  solve 
the  two  simultaneously.  This  is  done  in  the  obvious  manner:  at  each  point 
(i,j,k)  we  first  perform  the  divergence  relaxation  (3.9),  thereby  changing  the 
field  at  neighboring  points,  and  then  perform  the  Helmholtz  relaxation  (3.6), 
changing  the  field  at  (i,j,k). 

One  final  point  we  need  to  address  before  proceeding  to  examples  is  the 
question  whether  this  procedure  converges.   It  does  not.  Considering  just  the 
Helmholtz  equation,  even  when  all  the  boundary  values  are  specified  so  that 
there  is  a  unique  solution,  SOR  can  be  shown  not  to  converge  [24,25].   The 
reason  can  be  seen  heuristically  from  eqs  (3.5,3.6).   For  simplicity  consider 
Qpj  =  1  (Gauss-Seidel) .   Then  if  the  sign  of  4tt2A2  were  positive  in  the 
denominator,  SOR  would  replace  the  value  of  the  field  at  a  point  by  a  number 
somewhat  smaller  than  the  average  of  the  field  at  neighboring  points,  and  the 
process  would  converge.  With  the  negative  sign,  however,  the  field  at  the 
point  is  replaced  by  a  value  larger  than  the  average  of  neighboring  values, 
and  the  process  does  not  converge — the  fields  just  keep  growing. 
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Why  then   do  we  use   it?     There  are  two  reasons.      The   first   is   that   it   is 
still   possible  to  obtain  approximate  solutions   by  truncating  the  procedure  at 
some  point   even   if    it   does   not   converge.      If  we  monitor   the   volume  average  of 
the  residual,   we   can  use   as  our   result   the  field   configuration   corresponding 
to  the  minimum  average  residual.      This   truncated  SOR  will   be  used   in  the  next 
subsections,   with  some  success.      The  other  reason  for  using  this  rather 
suspect   procedure  is  that  there  is  not  an  obviously  better  one  available  which 
suits  our   needs — i.e.    does   not  require  inversion  of   a  singular  matrix   in  the 
problems  of  interest.      The    (scalar   or  vector)  Helmholtz  equation  just  happens 
to   be   awkward  numerically,    due  to  the  fact   that  when  it  is  written  in  the  form 
Ax   =  b    (cf.   the   discussion  preceding  eq   (3-5))    the  matrix  A   is  not    positive 
definite,    even  when  full   information   is   given  on  the  boundaries.      How  to 
overcome  this  difficulty  is  the  subject  of  current  research    (see  e.g.    [26]  and 
references   therein),    and   in  the  summary  we  shall   comment   upon   possible  future 
improvements   on  our  scheme.      For   the   present   examples  we  use  the   truncated  SOR 
method. 

3.2     Two-Dimensional  Examples 

In  this  subsection  we  present   two-dimensional   examples  which  were  used  to 

develop  the   procedure   and  test   dependence  on  various   parameters.      All   the 

examples   assume  a  cavity  of   rectangular   cross   section  with  perfectly 

conducting  walls  with  dimensions   axbxc.      The   TE£mn  mode  of   the   cavity  is 

given   by 

2,ttx        .      rrmy        .      nuz 

A     =  cos  sin  — ~     sin  , 

x  a  b  c 

S,b        .      &7TX  m-ny        .      nuz  ,  0    ir>\ 

A     = sin  cos  — r~     sin  ,  (3.10) 

y  ma  a  b  c 

Az   =  0, 
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where  we  have  chosen  unit  amplitude  and  assumed  m  *  0.   The  dimensions  will 
always  be  chosen  to  correspond  to  resonance  of  the  mode  used  to  specify 
"measurement"  results,  (t/a)2   +  (m/b)2  +  (n/c)2  =  4.   This  is  necessary  in 
order  to  have  any  solution  at  all  in  the  enclosed  cavity  in  the  absence  of  any 
source.   In  these  preliminary  examples  we  restrict  ourselves  to  cases  where 
only  one  component  (Ax)  of  the  vector  potential  is  nonzero,  and  then  in  order 
to  restrict  the  problem  to  two  dimensions  we  impose  constancy  of  Ax  in  the  x 
direction.  The  divergence  condition  (3.2a)  is  then  automatically  satisfied, 
and  we  need  only  treat  the  Helmholtz  equation  (3.2b). 

In  the  first  computation  we  assume  the  field  is  measured  at  one  point, 
the  center  of  the  cavity,  where  Ax  =  1.0,  Ay  =  Az  =  0.   The  cavity  dimensions 
are  taken  to  be  a  =  2b  =  2c  =  /2~7  which  corresponds  to  the  resonant  frequency 
for  the  TE011  mode.   The  tangential  components  of  A  are  fixed  at  zero  ateach 
wall,  and  Ax  is  required  to  be  constant  in  x.   The  relaxation  method  outlined 
above  is  used,  and  the  iteration  process  is  terminated  when  the  magnitude  of 
the  Ax  residual  averaged  over  the  volume,  <|RX|>,  reaches  a  minimum.   This 
first  example  is  atypical  in  that  <|RX|>  decreases  asymptotically  to  a 
constant,  and  so  our  solution  is  the  field  configuration  any  time  after  <|RX|> 
has  reached  its  limit  value,  which  is  about  0.4  x  10"1*.   A  value  of  Rjj  =  1.75 
resulted  in  the  fastest  convergence  to  the  limit. 

The  y  dependence  of  Ax  for  two  different  values  of  z  and  two  different 
grid  sizes  is  shown  in  fig.  3-1 .   The  solid  and  dashed  curves  represent  the 
field  configuration  for  the  TE011  mode,  which  is  the  only  mode  resonant  at 
this  frequency  which  is  consistent  with  the  measurement.   It  is  evident  that 
the  16x8x8  grid  is  sufficiently  fine  for  our  purposes.   The  z  dependence  of  Ax 
for  fixed  y  is  virtually  identical  to  the  y  dependence  for  fixed  z  and 
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therefore  is  not  shown.   The  final  results  for  Ay  and  Az  are  zero  everywhere, 
as  they  should  be. 

At  the  beginning  of  the  computation,  initial  values  of  the  field  at  all 
points  must  be  assigned.   Two  different  starting  configurations  were  tried. 
In  one,  A  is  initially  zero  at  all  points  except  measurement  points, where  of 
course  it  assumes  the  measured  values.   In  the  other,  a  smooth  initial 
configuration  is  generated  by  an  algorithm  which  interpolates  between 
measurement  points  and  boundary  points.   For  this  example,  the  same  results 
are  obtained  for  either  initial  configuration.   In  general  the  smooth  initial 
configuration  tends  to  yield  a  lower  minimum  <|RX|>,  anc^  we  snaH  use  it  in 
the  following  examples  unless  otherwise  noted.  We  expect  that  when  multiple 
solutions  exist  the  smooth  starting  configuration  will  result  in  the  smoothest 
of  the  possible  solutions.  This  should  correspond  to  the  lowest  mode. 

If  the  measurement  point  is  moved  from  the  center  in  this  TE011  example, 
the  quality  of  the  solution  deteriorates.   For  example,  if  the  measurement  is 
made  at  (y,z)  =  (b/4,  c/2)  then  the  computed  solution  at  the  center  is  about 
6%   higher  than  the  correct  solution,  as  shown  in  fig.  3~2. 

Turning  to  a  more  complex  field  configuration,  we  consider  the  case  of 

four  measurements  with  Ax  =  1.0  at  (b/4,  c/4)  and  (3b/4,  3c/4)  and  Ax  =  -1.0 

at  (3b/4,  c/4)  and  (b/4,  3c/4).   The  lowest  mode  consistent  with  those  values 

is  the  TE022  mode,  and  we  choose  a=b=c=  /^so  that  it  is  resonant.   In  this 

case  we  encounter  what  proves  to  be  the  typical  behavior  of  <|RX|>.   It 

decreases  with  successive  iterations  at  first,  but  it  levels  off,  begins  to 

increase,  and  eventually  blows  up.   The  value  ti^   =  1.75  again  results  in  the 

quickest  minimum,  and  other  values  of  Q#   do  not  result  in  a  significantly 

better  solution.  We  use  as  our  approximate  solution  the  field  configuration 

for  which  <|RX|>  is  a  minimum  with  fl^  =  1.75.   The  results  of  figure  3~3  then 
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follow.  The  solid  and  dashed  curves  represent  the  TE022  mode,  and  again  the 
agreement  is  very  good.  The  same  is  true  for  the  z  dependence,  which  is  not 
shown. 

Equally  good  results  were  obtained  using  four  measurements   representing 
the  TE021   mode,    using  a   cavity  with  dimensions  a=b=2c   =  /21.      The  measurement 
points   were  taken  to   be  the  same  as   in  the  TE022   case  above,   with  the  measured 
values   different  of   course.      Figure   3-iJa   compares   the   computed  Ax  to  the   TE021 
form   as   a  function  of   y  at    z=c/2,    and   fig.    3-^b  does   the  same  as   a   function  of 
z  for   y=b/4.      A  16   x    16    x   8  grid  was   used,   with  Jfy  =   1.75 

The  results   thus  far   have   been   very  encouraging.      Instructive 
difficulties   arise  when   a=b=c=  /0 . 8 ,'   however.      If  we  use  as  measurements   the 
TE021    field  at   two   points,    (y,z)    =    (b/4,    c/2)    and    (3b/4,    c/2) ,    the  agreement 
of   the   computed  results  with  the   TE021    field  shape   is  not   satisfactory.      This 
is   demonstrated  in  fig.    3~5.      In  this  and  succeeding  examples,   the  largest 
dimension   is  divided   into    16  divisions,   with  other   dimensions   divided 
proportionately.      Thus,    this   grid   is    16x16x16.      Away  from  the  measurement 
points   the   computed  field  falls  off   too  rapidly.      The   problem   persists  when 
the  number   of  measurement   points   is   increased  to  four,    as   is   evident  from  fig. 
3-6.      The   value  of   the  minimum  <|RX|>    ^s   somewhat   larger   than  the   TE022    and 
rectangular   TE021    cases:    0.25   x    10~2   as  opposed  to   0.20   x    10~2   and   0.11    x    10~2 
respectively. 

It   is  not   completely  clear  what   cause  underlies   these   difficulties.      The 

numerical  method  being  used  does   not   actually  converge  to  the  true  solution; 

it   approaches   it  and  then  wanders   off,   with  the   approximate  solution   being 

obtained  by  taking  the  fields  when   "closest"    to  the  true  solution    (as 

determined  by   <|RX|>)«      In  the   present   case    (TE021   with  square  cavity  cross 

section),    the   iterative   process  just   does   not   get   that   close  to  the   correct 
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result.   The  most  likely  reason  appears  to  be  that  there  are  two  possible 
resonant  modes  with  Ay  =  Az  =  0  at  this  frequency,  TE021  and  TE012,  and  the 
computation  does  not  pick  out  one  pure  mode. 

3.3    Three-Dimensional  Examples 

We  now  allow  variation  in  the  x  direction  and  treat  the  fully  three- 
dimensional  case.   This  requires  simultaneous  solution  of  both  the  vector 
Helmholtz  equation  (3.3b)  and  the  divergence  condition  (3.3a).   Again  we  start 
with  the  simplest  possible  field  configuration,  assigning  values  corresponding 
to  the  TE011  mode  at  measurement  points.   Six  measurement  points  were  used, 
located  at  (a/8,  b/8,  7c/8) ,  (a/4,  7b/8,  c/4)  ,  (3a/8,  b/2,  c/2) ,  (5a/8,  3b/4, 
3c/4),  (3a/4,  b/4,  c/4),  and  (7a/8,  7b/8,  7c/8).   The  points  are  scattered 
irregularly  to  cover  the  entire  volume  and  to  avoid  having  them  all  at  zeros 
of  some  low  order  mode.   The  box  dimensions  are  taken  to  be  2a=b=c=  1  //2]    and 
an  8x16x16  grid  is  used. 

The  computation  proceeds  in  much  the  same  manner  as  the  preceding  ones, 

except  that  now  the  grid  has  three  instead  of  two  dimensions,  and  the 

divergence  condition  relaxation  (3.9)  is  also  performed  at  each  point.   The 

quality  of  the  solution  depends  on  the  choice  of  the  relaxation  parameters,  % 

and  fty.   The  minimum  value  of  <|Rh|>  decreases  as  ft^  increases,  but  very 

slowly,  whereas  the  minimum  value  of  <|Ry|>  is  smallest  for  small  Q^   and 

increases  rather  rapidly  as  %  increases.   The  minima  of  both  residuals  are 

smallest  when  Qy  is  large.   The  optimal  choices  seem  to  be  around  Q^   =  0.2, 

fty  =  1.75,  which  are  the  values  we  used.   Representative  results  of  the 

computation  are  shown  in  fig.  3~7.   The  y  dependence  is  satisfactory,  but  only 

satisfactory,  and  a  similar  comment  applies  to  the  z  dependence  (not  shown). 

For  the  most  part,  the  x  dependence  is  very  good — the  computed  Ax  is  very 

28 


nearly  flat.  The  one  exception  occurs  along  the  line  through  the  center  of 
the  yz  plane.  Significant  deviations  from  constancy  occur,  particularly  at 
the  measurement  point  x/a  =  3/8  where  Ax  is  constrained  to  be  1 . 

The  TE011  results  are  not  bad  in  themselves,  but  they  do  not  bode  well 
for  more  complicated  cases,  which  will  have  greater  spatial  variation  and  more 
than  one  nonzero  component  of  A.   The  forebodings  are  borne  out  when  we 
usemeasurement  results  corresponding  to  the  TE022  mode.  With  six  measurement 
points  we  were  unable  to  obtain  acceptable  results.   In  fact  for  some  values 
of  £2fj  and  Qy   the  average  residuals  <|Rh|>  and  <|Ry|>  increased  steadily 
from  the  start  of  the  computation — the  best  field  configuration  was  the 
initial  guess.   Increasing  the  number  of  measurement  points  to  eleven  did  not 
improve  matters  significantly.   Also,  the  difficulties  persisted  when  we  used 
different  indicators  of  solution  quality,  such  as  max(R)  or  <R> . 

We  conclude  that  in  general  three-dimensional  problems,  given  just  a  few 
measurement  points,  solving  the  Helmholtz  equation  and  divergence  condition 
using  the  truncated  SOR  method  fails  to  obtain  useful  approximate  solutions. 
To  investigate  whether  the  numerical  method  is  more  successful  when  the  field 
is  specified  at  a  large  number  of  points,  we  did  computations  with  all 
components  of  the  field  fixed  at  all  boundary  points.   The  boundary  points  are 
chosen  because  they  can  be  thought  of  as  driving  the  solution;  when  the  field 
is  known  at  all  boundary  points  there  is  a  unique  solution.  We  wish  to 
investigate  whether  our  truncated  SOR  method  can  obtain  something  like  the 
correct  solution  and,  if  so,  whether  it  can  do  so  when  some  of  the  boundary 
points  are  not  given. 

We  have  obtained  results  for  a  number  of  low  order  modes.   A  good  example 

is  the  TE123  mode,  which  was  the  highest  mode  considered.   For  convenience,  a 

composite  residual  was  defined  by 
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R(a)  =   a  <|RR|>  +  (1  -  a)<|Ry|>,  (3.11) 

and  the  solution  was  chosen  as  the  point  at  which  R(a)  was  a  minimum.   A 
number  of  values  of  a  were  used,  and  the  results  were  relatively  insensitive 
to  the  exact  value.   In  the  computations  we  used  a  =  0.5.   Figure  3~8  compares 
the  computed  Ax  and  Ay  (Az  is  zero)  to  the  correct  TE123  results  as  functions 
of  y  at  x  =  a/4  and  z  =  c/8  and  c/2.   The  agreement  for  Ax  is  seen  to  be  very 
good  in  both  cases,  but  the  Ay  results  fail  to  reproduce  the  extremum  at  the 
center.   In  order  to  check  that  this  is  not  due  to  some  programming  clumsiness 
which  introduced  a  spurious  asymmetry  in  the  treatment  of  Ax  and  Ay,  we  also 
computed  the  TE213  mode  with  dimensions  a  and  b  interchanged.  We  found  that 
the  Ax  and  Ay  results  were  also  just  interchanged,  as  should  be  the  case.  The 
difficulties  apparent  in  fig.  3~8  can  therefore  be  taken  as  indications  of 
deficiencies  in  the  truncated  SOR  method  rather  than  the  consequence  of 
(obvious)  programming  blunders.   The  poor  results  at  central  y  will  also  be 
reflected  in  the  x  and  z  dependence  at  fixed  y.   This  is  seen  in  figure  3~9, 
where  the  (x,y)  =  (a/4,  b/8)  results  are  very  good  but  the  (a/2,  b/2)  result 
is  less  than  half  the  correct  answer. 

The  TE123  mode  constitutes  a  rather  demanding  test.   There  is  different 
spatial  variation  in  each  of  the  three  directions,  two  components  are  nonzero, 
and  the  components  have  different  spatial  dependence.   The  numerical  solutions 
are  quite  good  over  most  (-7/8)  of  the  volume  but  fail  in  the  central  region. 
A  finer  grid  led  to  only  slight  improvement.  When  the  fields  are  completely 
specified  on  the  boundary,  this  truncated  SOR  method  could  be  useful  in 
generating  a  starting  configuration  for  some  other  method  which  converges  (as 
it  was  used  in  [26]),  but  it  does  not  produce  satisfactory  results  by  itself. 
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Finally,    if   the   fields   are  specified  only  at   one-fourth  of   the   boundary 
points    (tangential   fields   still   vanishing  at   all   boundary  points),    no  minimum 
of  R   is  obtained  for   the  TE123   mode.      From   the  start  of   the   computation   each 
iteration  increases  R.      In  principle  we   could  fix  unmeasured  boundary   points 
by   interpolation  and  might   then  obtain  a  minimum  for  R,    but   at   this   time  that 
appears   to  be  beating  a  dead  horse. 

4.         DISCUSSION   AND   CONCLUSIONS 

We  have  suggested  an  approach  to  the   characterization  of  EM  environments 
which  is   based  on  the  numerical  solution  of  Maxwell's  equations  subject   to  the 
constraints   imposed  by   boundary  conditions   and  results   of  measurements   of  the 
field  at  a  relatively  small  number  of  points.      Some  success  was  achieved  in 
the  simpler  examples  used  to  demonstrate  the  approach,   but  results   in  more 
complicated  three-dimensional  examples  were  less  impressive.      The   comparisons 
made  were  of  the   field  configurations,   which   constitute  a  far  more   demanding 
test  than  does  a  global   quantity  such  as   the  average  power   density.      The 
results   of  a  less  stringent  test  could  well  be  less  disappointing. 
Nevertheless,    it   is   clear   that   the  method   is   not   yet  ready  for  real 
applications.      The  basic  obstacle  to  practical  use  of  the  approach   is   the   lack 
of   a  fully  successful  method  for  obtaining  numerical   solutions  of  Maxwell's 
equations   under  these  conditions. 

Two   different   numerical   techniques  were   investigated.      The  one  based  on 

the  action  functional   and  Hamilton's   principle  founders   due  to   our   inability 

to  numerically  locate  a  stationary  point   which   is  not   an   extremum.      Having 

previously  tried  a  Gauss-Seidel   type  of   procedure   [21,22]   and  having  found   it 

lacking,   we   performed  an   analytic  continuation  to   imaginary  time.      This   has 

the  advantage  that   the  action  then   has   a  minimum.      Furthermore,    the  solution 
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obtained  is  that  with  the  minimum  energy  density  consistent  with  the 
measurements.  This  technique,  however,  proved  to  be  inapplicable  to  cases  in 
which  there  was  any  spatial  variation  of  the  phase  or  any  phase  difference 
between  different  components  of  the  field.  We  therefore  abandoned  the  attempt 
to  develop  a  numerical  technique  based  on  the  stationarity  of  the  action.   If 
one  did  insist  on  pursuing  this  line  further,  one  possibility  would  be  to 
construct  and  numerically  minimize  the  function  defined  by 

2  2  2 

l        U3Av(i,j,k)J     l9A  (i,j,k)J     l3A  (i,j,k)J  J  *      ^*1J 
i,j,k  Y 

D2   is  a  nonnegative  function  whose  zeros   correspond  to  stationary  points  of 
the  action  and  therefore  to  solutions  of  Maxwell's  equations.      Obvious 
problems  with  that  tactic  are  that   the  construction  of  D2   is  quite  awkward  and 
that  there   is  no  guarantee  that   it  can  be  easily  minimized. 

The  second  technique  we  studied  in  detail  was  use  of  successive 
over-relaxation   (SOR)   to  simultaneously  solve  the  vector  Helmholtz  equation 
and  the  divergence  condition,   which  are  equivalent  to  Maxwell's  equations  in 
the   cases  considered   (sourceless,   free  space,   one  frequency).      Since  the  SOR 
method  does   not   converge  for   the  Helmholtz  equation,    the  sequence  of 
iterations  was   truncated  when  the  volume  average  of  the  magnitude  of  the 
residual  was  a  minimum.     This   produced  good  results  in  a  number  of   cases     but 
failed   in  others.      Even  when  the  field  was   completely  specified  on  the 
boundary,   the  results  were  not   always  satisfactory. 

There  are  various  other  possible  numerical  procedures  which  come  to  mind. 
Mittra  and  coworkers  at  Illinois  have  had  some  success  using  a  modified  Gauss- 
Seidel  method  on  problems  of  the  form  Ax  =  b  where  A  is  not  positive  definite. 
The  modification  they  make  is  to  introduce  a  convergence   factor  which  prevents 
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the   Gauss-Seidel   iterations   from   diverging    [27].      This   is   similar   to  our 
truncated  SOR  method,   with  the  truncation  replaced  by  a  gradual   damping   (and 
with  ft  =   1).      It  does  not   appear   likely  that   the  gradual   damping  would  cure 
the   ills  of  the  truncation  method.      The   divergence  problem   is  obviated   by  the 
truncation  method  too.      The   difficulty  is   that  the   iterative  procedure  does 
not  approach  very  close  to  the  correct  solution  in  some   cases,    and  the 
introduction  of  a  convergence   factor   does  not  address   that   problem. 

A  more  promising  alternative  would  be  to  try  a  truncated  conjugate 
gradient   procedure,   as  has   been  used  in  ill-conditioned  deconvolution  problems 
[28].      Some  work  would  be  required  to  adapt  the  procedure  to  our   case,    in 
particular  to  develop  a  criterion  for  when  to  truncate  the   iterations.      If  the 
method  can  be  applied,    it  could  prove  very  useful   since  proponents  of  the 
procedure  claim  it  yields  a  minimum-norm  solution.      For  our   problem  that 
corresponds   to  the  minimum  energy  density  solution. 

One  seemingly  relevant  method  which  proves  inapplicable  is  the  maximum 
entropy  method    (MEM),    as   reviewed  for   example   in   [29].      It  is   highly  effective 
for  reconstructing  images  which  have  been  blurred,    contaminated  by  noise,   or 
from  which  some  fraction  of  the   information  is  missing.      Our   problem  falls 
into  the  missing  information   category.      Unfortunately,    for  the  MEM  to  work   it 
requires  several   pieces  of   information  per  unit  cell,   where  the  cell   size 
corresponds  to  the  volume  over  which  the  information  at  a  point  has   been 
smeared.      Since  there   is  no   blurring  in  our   problem  there   is  no   information 
except   at   the  measurement   points,    and  the  MEM  is   powerless   to   fill   in  the 
empty  space. 

A  final   possibility  we  mention   is   to   put   the  Helmholtz  equation   in  the 

form  Ax   =  b,    as   described  in  Section   3,    and  then  multiply   by  the  adjoint  of  A 

(At)   to  obtain 
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Af  Ax   =   A+b.  (4.2) 

If   all   boundary  values   are  specified  the  SOR  method  will    converge  for  this 
equation  since  the  matrix  A  A   is   then   positive   definite.      Even  for 
incomplete   information  A  A  will   not   have  negative   eigenvalues    (although  some 
will  be   zero),    and  so  SOR  should  be  less   unstable.      As  with  D2  of  eq   (4.1), 
however,   A  A   is  not  so  simple  to  construct.      Furthermore,   the  rate  of  convergence 
for  iterative  solution  of  eq   (4.2)  would  be  expected  to  be  very  slow.      (This 
could  be   improved  by  using  SOR   on  Ax   =  b   to   precondition  x.) 

In  summary,   neither  of  the  numerical  methods  we  investigated  was 
successful   in  all   cases.      There  do   exist  other  methods  which  may  warrant 
further  investigation.     The  general  approach   (constrained  solution  of 
Maxwell's  equations)   sounds  attractive  for  characterization  of   complex 
environments,    but  its   practical   implementation  awaits  development  of  better 
numerical  methods. 
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Sources 


Fields  =  0 
on  surface 


Figure   2-1      Schematic  division  of  total   volume  into  volume  of   interest   (Vj) , 

buffer  volume   (VB) ,   and  remaining  volume  containing  all  sources    (Vs) 
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Figure  2-2  Rectangular  waveguide  dimensions  and  axes.  Also  depicted  is  the 
boundary  of  the  volume  over  which  the  action  integral  extends. 
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Figure  3-7     c)   Comparison  of  computed  results  to  TE0ll    field  for 

three-dimensional   computation.     Six  measurement   points  were 
used, As  a  function  of  x  for  various   y  and  z. 
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